CONSTRUCTING GENERALIZED SYNCHRONIZATION 
MANIFOLDS BY MANIFOLD EQUATION 
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Abstract. Full understanding of synchronous behavior in coupled dynamical systems beyond the 

identical case requires an explicit construction of the generalized synchronization manifold, whether 

... we wish to compare the systems, or to understand their stability. Nonetheless, while synchronization 

^\ ' has become an extremely popular topic, the bulk of the research in this area has been focused on the 

C " 3 identical case, specifically because its invariant manifold is simply the identity function, and there 

f^^ , have yet to be any generally workable methods to compute the generalized synchronization manifolds 

^\] ■ for non-identical systems. Here, we derive time dependent PDEs whose stationary solution mirrors 

' exactly the generalized synchronization manifold, respecting its stability. We introduce a novel 

1-" ' method for dealing with subtle issues with boundary conditions in the numerical scheme to solve the 

yJ ' PDE, and we develop first order expansions close to the identical case. We give several examples 

^Jh , of increasing sophistication, including coupled non-identical Van der Pol oscillators. By using the 



in 



Q 



^ 



(N 



>< 



manifold equation, we also discuss the design of coupling to achieve desired synchronization. 

Key words, generalized synchronization, computational invariant manifold, coupled dynamical 
systems 



^^ i 1. Introduction. Following the fundamental finding in [2] that even chaotic 



systems can synchronize through coupling, synchronization phenomena in coupled 
oscillator systems have been studied extensively during the past few years, in the con- 
text of dynamical systems ([71[I71[l[l[Sl[TSl[51[Hl[Ii]), control theory {[TD]), complex 
networks ([TTl [191 HO]), laser physics ([15l [18]), etc. Besides identical synchroniza- 
tion, studies have been carried out to describe other types of synchronization, such as 
h> ' generalized synchronization O HI [6l [18] , phase synchronization O [8] , and anticipated 

f^ I and lag synchronization [H [15] . 

OO . The simplest form of synchronization between two oscillators is the complete 

(or identical) synchronization, meaning that after a transient time, the difference 
between the two oscillators approaches zero. However, a more generalized form of 

'bj I synchronization can often occur between two non-identical systems. 

In the combined state space of such systems, an invariant manifold consisting 

f~^ I of the synchronized states, termed generalized synchronization manifold, is a central 

object of study in this context, and its form ([H]) and stability ([9l [16]) are two 
important properties besides its existence. The former tells the relationship between 
corresponding dynamic variables in the synchronization state, and the latter deter- 



JH I mines whether a trajectory near the manifold evolves closer or farther away. While a 

great deal of work in synchronization focuses on identical synchronization, primarily 
because the synchronization manifold in this case is described simply by the iden- 
tity function, much less work (mainly on the detection of existence and stability of 
generalized synchronization, see O [4] [6]) has been carried forth on the ubiquitous 
generalized synchronization phenomenon, specifically because there has been little 
work to explicitly construct the corresponding synchronization manifold. One possi- 
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bility, as reported in [13] , is to approximate the manifold using orthogonal functions 
expansion. 

Although one can simply run the system to find a portion of the generalized 
synchronization manifold that is covered by the attractor of the coupled system, the 
study of the entire manifold is still necessary and important. One reason is that the 
transformation from one state variable to the other that defines the manifold contains 
information on the correlation between the dynamics of the two systems in the whole 
state space, not just on the attractor. Furthermore, the manifold outside the attractor 
may embed relevant unstable invariant sets, which can become stable as parameters 
change. 

Our recent work strives to fill this gap. In this paper we present a PDE (manifold 
equation) which describes the invariant manifold of generalized synchronization be- 
tween two coupled oscillators as well as a variational equation regarding the transverse 
stability of the manifold. We also develop a time dependent PDE whose stationary 
solution is a generalized synchronization manifold, and present straightforward nu- 
merical schemes to solve the equation. Several examples, such as ID and 2D coupled 
oscillators and coupled Van der Pol oscillators are used to illustrate our method of 
constructing the invariant manifold under generalized synchronization. The design of 
coupling in coupled system to achieve desired synchronization is also discussed by the 
use of the manifold equation. 

In Section II, we give background regarding generalized synchronization, leading 
to the manifold equation together with discussion of stability of the manifold and 
perturbation analysis relative to the identical case. In Section III we show a specific 
form of the manifold equation regarding coupled oscillator system with almost the 
same individual dynamics and a time dependent PDE used to solve for the original 
equation. In Section IV we show some examples of invariant manifold. In Section V 
we show how to use the manifold equation to design coupling function between two 
oscillators to obtain desired form of synchronization. 

2. Equations for synchronization manifold. For the discussion in this pa- 
per, consider the following equations to describe the dynamics of two non-identical 
oscillators wi and W2'- 

wi = f{wi,W2,^J.l), 
(2.1) W2= g{wi,W2,y.2)- 

here wi e 3?™, W2 e 3?", / : 3?™ x 3?'" x 3?! ^ K™ and 5 : gj™ x 3?™ x 3?i ^ 3ff™ 
where f £ C^ and g G C^ (both / and g have continuous first order derivative) . Note 
that the coupling functions which enable the two oscillators to synchronize, have been 
included in the general form of the functions / and g. /ii and 112 are considered as 
parameters of the functions. 

In this paper, we will focus on perturbations from identical synchronization. In 
such a case, we require that f{wi, wi, ni) = g{wi, wi, /ui), meaning that when the two 
systems have the same parameters and same value of variables, they have the same 
individual dynamics, often with no effective coupling between them. 

2.1. Manifold Equation. Let (wis, W2s) denote the synchronization state. The 
synchronization manifold at synchronization state W2 — ^(wi) istheset {{'Wi,W2)\w2 = 
<i?(wi)} where $ is a time independent transformation. In this paper we will abbrevi- 
ate this notation and just use the relationship $ to indicate the form of this manifold. 
To find the synchronization manifold of this coupled system, we mean to find some 
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C^ function $ ($ : W -^ K™) such that 

(2.2) W2s = $(m.). 

Suppose such $ exists. 

Then, along the synchronization manifold, 

(2.3) u;2. = '^^^ = D^wu) ■ w\s 

at 

Using Eq. (|2.ip in Eq. (12. 3|) for the synchronization state, and replacing wis by 
wi for convenience, we obtain the following manifold equation: 

(2.4) D^{w) ■ fiw, <^iw),fii) = giw, $(w), /is). 

This is a partial differential equation with unknown function $ ~ ($i, ..., $„) which 
is of dimension m and each component $i has to variables, and has been used for the 
study of generalized synchronization pTSj and center manifold approximation [1.^. 

In the case that fJ-i — 1^2 = fJ- (coupled identical oscillators), a particular solution 
of Eq. (|2.3p is $(w) = w, since f{w,w,^) — g{w,w,fi) which describe the identical 
manifold under complete synchronization. However, in general cases, it is much harder 
to find a solution to this PDE. We will leave it to the following to discuss some 
appropriate approaches in obtaining the solution. 

2.2. Stability of the Manifold: Variational Equation. In addition to the 
form of the manifold, the stability plays another important role in the study of syn- 
chronization. To measure the stability of the synchronization manifold W2s = ^{wis) 
obtained by Eq. (|2.4[) . define 



(2.5) ^ = W2-^{wi) 

to be the transverse perturbation from the manifold. Then from Eq.s (|2.1[) . (|2.3[) and 
(1231), we have 



^ — W2 — D^{wi)wi 

= giwi,W2,tl2) - D^{wi) f {wi,W2, y.i) 

= g{wi, $(u;i) + C, M2) - .9(^1, $(wi), /i2) + g{wi,'^{wi),fi2) - -D$(wi)/(u;i, $(u;i), Mi) 
i'M^{wi)f{wi, $(u;i), fii) - D^wi)f{wu'i>{wi) + ^ Mi)- 



The term g{wi,^{wi), ^2) - D^iwi)f{wi,^{wi),fii) = by the manifold Eq. (1231) . 
so expanding the other terms to the first order in ^, we obtain the local variational 
equation: 

(2.7) i= [D^^g{wi,W2, tJ.2) - D<^{wi)D^J{wi,W2, fJ.i)]\^^=<s>(wi)C 

This variational equation can be used to study the local stability at the points on 
the invariant manifold. However, we note that when this equation does not give a 
uniformly asymptotically stable solution, the asymptotic stability of synchronization 
may not be correctly predicted from this equation, see examples given in [14]. From 
Eq. (|2.71) we again see the importance of the knowledge of the manifold W2s — ^{wis)- 
Although there are approaches discussing the stability of the synchronization man- 
ifold without knowing explicitly the form of the manifold, it will be more natural 
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and convenient to study the stability directly with the knowledge of the form of the 
synchronization manifold. 

In this paper, we will focus on the form of the invariant manifold. The detailed 
discussion regarding the stability of the invariant manifold under generalized synchro- 
nization will be reported in the future work. 

2.3. Specific Form of Manifold Equation for Perturbed Oscillators. 

Suppose that the individual dynamics for the two oscillators wi and W2 in isola- 
tion are almost the same, but with a small difference. To be explicit, we mean that 
in Ea. l2.1l there is a small mismatch in the parameters: fi2 — l^i + f ■ Then, one would 
expect the synchronization manifold to be close to the identical synchronization man- 
ifold and $ to approach the identity function as e — > 0. We empirically expect the 
synchronization manifold to have the form 

(2.8) W2. = $K.) = wu + eH{wis) + 0{e^). 

The main interest will be to find H and thus we will have a first order approximation 
for the synchronization manifold. Following the similar procedure to the derivation 
of the manifold Eq. (|2.4p . after neglecting higher order terms we obtain a PDE for H: 

(2.9) +D^^^[g{wl,W2,^J.l) - fiwl,W2,^J-l)]\w2=wlH{wl). 

3. Solving the Manifold Equation in Practice. In the previous section we 
derived the general equation for the synchronization manifold of non-identical oscil- 
lators as well as the corresponding variational equation for the study of stability of 
the manifold. We conclude that in order for such transformation W2 — $(wi) to exist 
when t ^ oo, (i.e., $ is a stationary solution of Eq. 13. 21 when t is sufficiently large) <& 
has to satisfy Eq. (12. 4|) . So the problem of finding the invariant manifold reduces to 
finding the solution of the PDE (12. 4p with appropriate boundary conditions. 

However, the boundary condition for Eq. (j2.4p is not easily attainable without 
actually solving the original ODE system. There may exist more than one synchro- 
nization manifolds in general, if none of them is globally stable. Our main interest in 
this paper is to use the perturbed manifold equation (|2.8p to find the closest one to 
the identical manifold, which tells us how the form of synchronization changes if the 
two oscillators are not perfectly the same. 

3.1. A Time Dependent PDE for the Manifold. In order to develop a 
workable numerical scheme to solve the manifold equation (|2.4p . we first derive an 
evolution equation for a mapping relating the coordinates wi and W2 of the two sys- 
tems obeying (|2.ip . The idea is that if we start with a reasonable initial guess for the 
mapping and evolve it according to the dynamics of the coupled system, the mapping 
would approach the synchronization manifold, provided that the manifold is asymp- 
totically stable. Indeed, the synchronization manifold turns out to be a stationary 
solution of the evolution equation, whose stability mirrors exactly the stability of the 
synchronized systems. 

We assume the existence of a smooth time-dependent mapping from the phase 
variable wi to W2, i.e., 

(3.1) W2 = (j){wi,t), 
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as an ansatz. Differentiating this with respect to t and using Eq. (|2.ip . we see that 
(t>{wi, t) satisfies 

(3.2) — + -^f{wi,(l),fii) ^g{wi,(j>,fi2), 

which is a PDE that describes the time evolution of the functional relationship be- 
tween the states of the two systems. If we set -^ = 0, the equation reduces to 
Eq. (|2.4p . so a synchronization manifold, if it exists, is a stationary solution. Since 
the evolution equation (|3.2[) reflects the co-evolution of the coupled systems, we ex- 
pect a synchronization manifold to be asymptotically stable with respect to Eq. (|3.2p . 
i.e., 

(3.3) $(wi) = lim (/)(wi,t), 

t — >oo 

whenever synchronization is asymptotically stable, and the initial condition for Eq. (j3.2p 
is within the basin of attraction. 

If we look for a synchronization manifold of form (j2.8p . then we can write the 
time-dependent manifold as 

(3.4) 0(wi,i) = wi +eh{wi,t) +0(e^) ~ wi +€h{wi,t). 

Neglecting higher order terms in e we can derive the corresponding time dependent 
PDE for the function h{wi, t): 



-^ ^D^,^_g{wi,wi, ^J.2)\t,2=p.i + -DtuJg(wi,u;2,Mi) - f{wi,W2, t^i)]\w2=wih 



dh 
'dt 

(3.5) --^f(wi,wi,fj.i). 

The stationary state of h gives H corresponding to the synchronization manifold: 

(3.6) H{wi) == lim h{wi,t). 

t — >Ctt 

3.2. An Iterative Scheme to Find the Stationary Manifold. One ap- 
proach to solve Eq. (|3.5p is by successive approximation: start with some initial form 
and evolve the solution according to the equation until the approximation gets close 
enough to the true solution. To obtain an iterative scheme to approximate the syn- 
chronization manifold, we first discretize t in (j3.4p by considering t = tQ,ti, ... where 
Iq — Q and t„=i — tn '■= Tn- Let us use /i„(u'i) to represent h{wi,tn). Also, for 
convenience, let 

(3.7) b{wi) -.^ Df,.^g{wi,wi, ^i2)\^L2=^n, 

(3.8) B{wi) := D^^[g{wi,W2,fJ.i) - f{wi,W2, ^J.l)]\^u2=wl■ 

Then the iteration scheme can be described as 

• Given initial guess /iq. 

• For n = 1, 2, ..., until convergence, do 



(3.9) /i„ = /i„-i +T„_i 



b{wi) + B{wi)hn-i TT — fiwi^wi,fii] 

OWi 
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To analyze the convergence of this scheme, let h* denote a fixed point of our 
iteration operator, i.e. 



(3.10) /i* = /i*+T„_i 

Then we have 



dh* 
b{wi) + B{wi)h* - - — /(wi,wi,/ii) 

OWi 



hn - h* = hn-i - h* 



Tn-1 



B{wi){hn-i - h*) - 



d(hn-l 



h*) 



(3.11) 


so that 


\\K-h*\ 


(3.12) 



{Im + Tn-lBw^){hn-l — h*) — T„_i 



<I|I„ 



Tn-lBwi\ 



dwi 

g(fe„_l - h*) 

dwi 



d{hn 



-f{wi,Wl,fll) 



f{wi,Wl,fll) 



h*) 



dwi 



||/(wi,wi,Aii)| 



for appropriate matrix and vector norms. 

The relevance of each term in this error propagation equation is interpreted as 
follows. The first term ||/„j + t„-iBw^\\ ■ \\hn-i — h*\\ has propagation factor ||/„j + 
T„_ii?u,j|| smaller than 1 if synchronization is stable, while the second term |t„_i| • 

^ "c)^~ — - ■ ||/(wi,wi)|| corresponds to the total variation in space of the error 
function, which in general is not guaranteed to converge to 0. 

The difficulty of controlling the variation part in the second term above suggests 
the following adaptive step size scheme: 



(3.13) 



K 



hn-l +Tn-1 



b{wi) + B{wi)hn^i -a„-i^7^^^ — f{wi,wi) 

OWi 



where a„_i S[0,1] is a control factor that controls the contribution of the total vari- 
ation to the next iteration: when the total variation is small, choose a„_i close to 
1, and when the total variation is large a small, Q!„_i will be preferred. A necessary 
(but not sufficient) condition for the solution to converge is that a„_i ^ 1 as n — > cxd. 

3.3. Spatial Discretization in 2D. To be specific, we will demonstrate the 
2D discretization in space. The ID case can be obtained easily from the 2D scheme 
and for higher dimensions our scheme can be modified to work. 

Suppose we have the 2D version of (|3.9p . and denote wi — {x,y)\ and h{wi) = 
{u{x,y),v{x,y)y. Also, let 



(3.14) 
and B('Wi) 



h{wi) = {bi{x,y),b2{x,y)y, 
f{wi,wi) = {fi{x,y),f2{x,y)Y 



Bii{x,y) Bi2{x,y) 
B2i{x,y) B22{x,y) 
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Then the iterations for u and v wiU be 



U„ = U„-l +T„_l 



(3.15^„ = u„_i + r„_i 



h + [BiiUn-l + Bi2Vr,-l) ~ I — /i H ^ ,/2 

02 + (-B2lW„_l + B22Vn-l) - I — ^ /l H ^ 72 



given initial guess Uq and t;o. Suppose we want to solve the equation in a domain 
[a, 6] X [c, d] , we can form a grid on this domain with mesh size A^, and Aj, respectively: 



(3.16) 



Xi = a + iAx, i = 0, 1, 2, ..., Nx, 
y, =c + jAy, j = 0,1,2,.. .,Ny, 



where N^ ~ ^^ and Ny = ^-^. By using the central difference operators, we have 

dUn-i{i,j) _ Un^i{i + l.j) -^„_l(^ - l,j) 
dx ^ 2A^ 

dUn-l{i,j) Un-lii,j + 1) - Un~l{i,j ~ 1) 



dy 



2A„ 



(3.17) 



dVn-l{i,j) _ «n-l(z + l,j) -Vn-l{i - l,j) 

dx "* 2A^ 

dVn-l{i,j) Vn-l{i,j + 1) -Vn-l{i,j - 1) 



dy 



2A„ 



the iteration scheme for 2D problem (|3.15p at each grid point {xi, yj) will be 



Un = Un-1 +T„_i 



bl + (BllU„_i + Bi2Vn-l) 



Un-l{i+l,j) - Un-l{i-l,j) Un-l{i,j + 1) -Un-l{i,j -I) 

/I H ;t7 /2 



Wn = «„_! + r„_i 



2A, 

&2 + {B2lUn-i + B22Vn-l) 



2A„ 



,g -^gN I' vu^iji+^J) -Vn-iji- l,j) , «n-i(i,j + 1) ~««-i(i,j - 1) „ 



2A, 



2A„ 



Here the functions without arguments are evaluated at grid points (xi,yj). 

3.4. On Boundary Condition and Computation. There is one issue for 
implementing the iteration p.lSp : the values on the boundary of the domain is difficult 
to compute in general. One approach to avoid this issue is to use dynamically shrinking 
domains. Start by assigning values at all points in a initial domain that is much 
larger than the domain of interest. Then the iteration is defined properly at all 
interior points, but not on the boundary, so wc make the domain smaller by discarding 
all boundary points for the next iteration. Repeating this process, we produce a 
sequence of domains of decreasing size on which successively better approximation of 
the synchronization manifold is obtained. By taking the initial domain to be large 
enough so that the iteration converges before the domain becomes smaller than the 
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original domain of interest, we obtain an approximate solution on the entire domain 
of interest. 

Intuitively, this method works by choosing the boundary of the computational 
domain to dynamically change with the propagation wave front of the effect of the 
boundary of the initial domain. Since the solution at a given fixed point in the domain 
will not be affected by the choice of values on the boundary of the initial domain until 
this wave front reaches the point, any point that remains in our shrinking domain 
should not be affected by the choice boundary condition on the initial domain. The 
influence of the unknown boundary condition of the outer domain becomes progres- 
sively less important as time progresses and the inner domain shrinks. 

4. Examples of Invariant Manifold. As we discussed in the previous sec- 
tions, the manifold equation (|2.4p is important because it allows further study of the 
synchronization between coupled oscillators. In this section we show several examples 
of constructing a generalized synchronization manifold. In particular, we consider 
examples in which the dynamics of the two oscillators are different but only with 
minor difference. A motivation for the study of this special case is that in any real 
dynamical systems no oscillator is physically perfect in the sense that the parameters 
or even the structure we use to mold these oscillators bear some small error. 

4.1. A Coupled ID Oscillators. One-dimensional examples provide an ex- 
cellent starting point of investigation for studying invariant manifolds for higher- 
dimensional systems. Consider a simple case where 

(4.1) i=l, xeSR. 

We couple this oscillator with another one, but with some perturbation function on 
it so that the two oscillators are not identical: 

i — 1, 

(4.2) y^l + {x-y)+e{p{x)+p'{x)). 

where p is any smooth function and p' denotes its derivative. Here the (x — y) term 
represents the one-way coupling from x to y and the e{p{x) +p'{x)) term represents 
some small perturbation on the y oscillator. Note that when e = the synchronization 
manifold is y — x, which is stable. This oscillator is designed to have <I>(x) — p{x) 
as a stable invariant manifold, which can be checked by introducing an appropriate 
change of coordinates. This type of oscillator always results in an intrinsic relationship 
between x and y subcomponent, given by y = P^x). 
For example, if let p{x) = sin(x), then we have 

i;= 1, 

(4.3) y = 1 + {x — y) + e(sin(a;) -|- cos(a;)). 

Now for e ^ 0, the invariant manifold <I>(a;) is determined by the manifold equation 
(1231), which yields 

d^(x) 

(4.4) — ^ = -$(x) + {l + x) + e(sin(.T) + cos(.t)). 

ax 

This equation is a simple ODE that can be solved analytically using the multiplication 
factor method. The solution is 

(4.5) $(a;) ^ x + e sin(2;) + Ce~^ 
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where C is a constant that is determined by the boundary condition. In this family of 
invariant manifolds, the only one that corresponds to the perturbation of the identical 
synchronization (<&(x) — > x when e — > 0) is obtained when C = 0: 

(4.6) $(x) = x + esin(a;). 

The stability of the manifold can be determined by (|2.7p , in this case we obtain 

(4.7) e = -e, 

so for t —>■ cxD, the invariant manifold is exponentially stable. 

The manifold can also be obtained by using the iteration scheme we described in 
the previous section: suppose the form of the manifold ^(x) — x + eH{x). To find 
H we look for the corresponding stationary solution of the PDE described by (|3.5p . 
Following the iteration scheme p.9p . we start with /ig = and choose r„ = r = 0.1, 
and obtain the following sequence of functions: 

hi ~ 0.1 • sin(a;) + 0.1 • cos(a;), 
/i2 = 0.2 • sin(a;) + 0.18 • cos(a;), 
/iio = 0.8340089600 • sin(a;) + 0.3315041568 • cos(a;), 

/iioo = 0.9999965626 • sin(x) - 0.00004893563553- cos(2;). 

As we can see, this sequence of iteration functions converge to sin(a:), which is the 
true form that appears in the manifold. We plot the first 20 iterates in Fig. 14.11 to 
visualize the convergence of iterations. 

4.2. Coupled 2D Oscillators . We now consider two coupled oscillators and 
study the invariant manifold that describes the long-term relationship between them. 
More specifically, consider a unidirectionally coupled system made up of two different 
oscillators: 

xi = 1, 

yi = 1 + [xi -yi) +ei{pi{xi) +pi(a;i)), 

X2 = 1 + k{xi - X2), 

(4.8) y2 = 1 + {X2 ~y2) +I^2{P2{X2) +P2(a;2))- 

Here the first oscillator is sending signal to the second through the x component, 
where k is the coupling strength. For large enough k we expect the corresponding 
X components to synchronize identically, while the form of synchronization of the 
y component can depend on the choice of e and p. However, note that when x 
synchronizes, the coupling term becomes small and we have two individual oscillators 
which we know how to solve. For t — > 00 we expect to have the following relations to 
hold, at least in the approximate sense: 

a;2 = a;i, 

2/1 = xi +eipi{xi), 

(4.9) V2=X2+e2P2{x2)- 
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Fig. 4.1. Successive approximation of the manifold form H(x) corresponding to the system 
i4-^ ! here ht(x) corresponds to successive functions obtained from (3. 9|) . voith step size r = 0.1. 
The approximation tends to the true form H{x) = sin(a::). 



Let e :— 62 — ci measure the order of mismatch of the two oscillators, and we look 
for the relationship y2 — yi + eH{xi,yi). Plugging Eq. (|4.9p into this equation, we 
obtain 

(4.10) eH{xi,yi) = £2^2(2:1) - eipi(a;i). 

Thus, we expect the invariant manifold between these two oscillators to be given by: 



(4.11) 



a;2 = xi, 

y2^yi+ <^2P2{xi) - £iPi{xi). 



Indeed, one can show that these equations do give an exponentially stable synchro- 
nization manifold by introducing an appropriate change of variables. 

As an example, let us consider pi{x) — P2{x) — sin(x), but with ei 7^ 62. Setting 
e = £2 — ei, the invariant manifold is given by 

X2 = Xi, 
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(4.12) 2/2 = yi + esin(a;i). 
To check the iteration scheme in this case, let 

X2 = xi +eu{xi,yi), 

(4.13) y2 = yi + ev{xi,yi) 

be the form of invariant manifold and apply the iteration scheme p.9p with h{xi, yi) = 
{u{xi,yi),v{xi,yi)y. 

Starting with the initial functions ua = vq = 0, wc find that 

uq = 0, t;o = 0, 

ui = 0, wi = 0.1 • sin(a;) + 0.1 • cos(x), 
U2 = 0, W2 = 0.2 • sin(a:) + 0.18 • cos(x), 
uio = 0, wio = 0.8340089600 • sin(a;) + 0.3315041568 • cos(a;), 

uioo = 0, fioo = 0.9999965626 • sin(a;) - 0.00004893563553 • cos(a;), 

and Un converges to while u„ converges to sin(xi), as expected. 

4.3. Coupled Non-Identical Van der Pol Oscillators. The equation of Van 
der Pol oscillator can be written as a first order system: 

x = y, 

(4.14) y = -X + fi{l - x'^)y 

where x represents position, y velocity, and fi is the strength of the nonlinear damping. 
For /i <C 1, there exists a stable limit cycle of radius 2. Now consider two of these 
oscillators, with slightly different values of the parameter /i, and coupled through the 
X component from the first to the second. The equations describing this system are: 

xi = yi , 

m ^ -xi + fi(l -xl)yi, 

X2 = 2/2 +k{xi - X2), 

(4.15) y2^-X2 + {fi + e){l~xl)y2. 

Here in the second oscillator the k(xi — X2) term represents the one-way coupling 
from the first oscillator, and e <C ^ is the mismatch in the damping strength. Since 
here we have two different individual systems coupled together, the identical manifold 
is not invariant under this dynamics. However, for large enough k, when i — > 00 we 
numerically observe that the difference of xi and X2 is O(e^) while the difference of 
yi and y2 is 0(e). In Fig. 14. 21 we show an example of the deviation of trajectory from 
the identical manifold, with parameters chosen to be /i = 0.1, e = 0.01 and k — 20. 
Generalized synchronization appears instead of complete synchronization in this case. 
To demonstrate the effectiveness of our numerical scheme, use the 2D iteration 
scheme (j3.15p . start with initial guess close to the theoretical prediction u{x,y) = 
0, v{x, y) = x—^x^ (see Appendix for derivation), and compute the iterations starting 
on the domain: [—5, 5] x [—5, 5] with mesh size 0.005 and time spacing r — 0.0005. 
We terminate when the domain shrinks down to [—2.5,2.5] x [—2.5,2.5] and take the 
solution at that time to be the approximated stationary solution. In Fig. 14.31 we 
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Coupled Van der Pol Eqs AX vs. time(linear-log plot) 




10" 



10" 



10" 



>- in-4 



250 300 350 400 

t 

Coupled Van der Pol Eqs AY vs. time(linear-log plot) 

I I I 



< 



10 



10" 



10" 



10 



200 



250 



300 
t 



350 



400 



Fig. 4.2. Deviation from the identical manifold. The time series comes from the system f4-15\ l 
starting from (a;i(0),3/i(0), X2(0), j/2(0)) = (1.5,1.5,1.5006,1.5107). Here in the two panels we plot 
the intervalt G [200,400]. The curve in the upper panel shows the difference AX{t) := |x'i(t) — a:2(i)| 
and in the lower panel we plot AY{t) := \yi{t) — y2{t)\- Red line in the upper panel shows the value 
e^ = lO"'* and the black line in the lower panel indicates e = 10~^. 
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plot the solutions u and v, with a typical trajectory obtained from the original ODE 
system. 

Figure [44] shows the successive error Cn = max(||u„ — Mn-i||oo, H^n — "yn-iHoo) 
measured by the infinite norm. We see exponential decrease in this error, which in- 
dicates that the iterative scheme is converging. In Fig. 14.51 we plot the distance of a 
typical trajectory to the manifold obtained by this scheme, which suggests that the 
computed manifold indeed gives a good first-order approximation for the synchroniza- 
tion manifold. 

5. Design Coupling to Satisfy Desired Synchronization. In many real ap- 
plications it is useful to have the knowledge of how to design specific coupled system 
so that the resulting behavior satisfy certain purposes. Mathematically we are seek- 
ing for a method to design the coupling between two general oscillators so that the 
resulting invariant manifold will have the prescribed form. In this section we discuss 
how to achieve the above goal by making use of the manifold equation we derived 
earlier. 

Assume that without coupling we have two separate oscillators described by the 
following equations: 

wi = F{wi), 

(5.1) W2 - G{W2). 

where wi G 3?™, W2 G 3?", i^ : SR" ^ SR" and G : K" ^ K™ where F & C^ and 
G G C^. Now suppose we want wi and W2 to synchronize, and furthermore, the 
synchronization manifold is prescribed as W2 — $(wi), the question is how to couple 
the two oscillators so that the resulting invariant manifold will have the desired form 
and be stable. 

5.1. Form of Coupling Function. We first solve for the coupling function so 
that the corresponding coupled system will have the invariant manifold satisfying the 
prescribed form. For convenience, we use a one-way coupling, so that the system 
becomes 

wi = F{wi), 

(5.2) W2^G{W2)+K{WI,W2). 

where K : K™ x 5R™ -^ 3?™ represents the coupling from x to y. In order not to make 
the question too broad, we assume that the coupling function has the form 

(5.3) Kiwi,W2)^a{L{wi)-H{w2)) 

where L : 3?™ -^ 3?™ and i? : 3?™ ^ 3?™ both G C^ . cr is a scalar which measures the 
strength of coupling, often called the coupling coefficient^ or coupling strength. Then 
in order that the synchronization manifold of the corresponding coupled system 

wi = F{wi), 

(5.4) W2 = G{w2)+a{L{wi)~H{w2)) 

has the desired form W2 — $(wi), by equation (j2.4p . we have: 

(5.5) D<^{wi)F{wi) = G($(wi)) + (t{L{wi) - H{<P{wi)). 
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Fig. 4.3. Numerical solution of the time dependent PDE corresponding to the system f^. J5| ). 
In the upper panel we plot u, and the black curve landing on u is obtained from plotting time series 
(xi{t),yi{t), i£2LJ — 2LLLi.\ obtained from original ODE. Similarly, in the bottom panel we plot v as 
well as the time series (xi{t),yi{t), ^2L2 — ELLlj. 
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Fig. 4.4. Successive error e-n, in the 2D iteration scheme for the time dependent PDE of the 
non-identical coupled Van der Pol system. Here the error is measured by the infinite norm. 



Note here that we have two functions L and H to choose while the manifold only 
gives one constraint, we need to fix one of the coupling forms. Suppose we are given 
the form of H{w2), and want to find L(wi). Solving for L{wi) yields 



(5.6) 



1 



L{wi) = -[£>$(u;i) • F{wi) - G($(u;i))] + i/($(wi)). 



So for any given form of the synchronization manifold W2 — ^(wi), we are able to 
design the corresponding coupling between the two oscillators. 

Notice that for identical oscillators (F — G), if we want ^{wi) — wi (identical 
synchronization), then 



(5.7) 



D^{wi) = 1^, 



is the m X m identity matrix so the first term in Eq. ()5.6|) disappears and the second 
term is just H{^{'Wi)) = H{wi), so we have 



L{w,)^H{wi), 



(5.8) 



meaning that for the coupled system (|5.4[) to have the identical synchronization man- 
ifold, the coupling functions L and H should have the same form, as expected. 

5.2. Adjusting Coupling Coefficient to Obtain Stability. In Eq. (|5.6p we 
have the form of the coupling function which guarantees that the coupled system 
described by Eq. (|5.4p have the invariant manifold W2 = $(wi). The correspond- 
ing variantional equation along the manifold can be obatined using Eq. (12. 7p . which 
becomes 



(5.9) 



i = D[Giw2) - (rH{w2)]\w2='S>(wi)C 
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Coupled Van der Pol: e (t) vs. t 
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Fig. 4.5. Distance of a typical trajectory to the manifold obtained by numerical scheme. Here 
ex{t) '■= \^2{t) — xi{t) — eu{xi(t), 3/i(i))| and ey(t) := \y2{t) — yi{t) — ev{xi{t), yi{t))\ where u and v 
are solutions by numerical scheme on grid points and for points in between grid points we use spline 
interpolation to obtain its corresponding value. Here we see that Cx ~ o(e) and Cy ~ o(e), suggesting 
that the first order approximation in e (e = O.OlJ is achieved. 
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Stability can be obtained by altering the value of a^ the coupling strength, if it is 
possible to have stable synchronization. Furthermore, with this variational equation 
as well as the form of the invariant manifold (known in advance), we are able to 
analyze the stability of not only the whole manifold but also individual points or 
subsets on the invariant manifold. 

5.3. Examples of Application. One case of the design of coupling form ap- 
pears in linear systems with perturbations. Suppose we have two uncoupled linear 
systems: 

w\ — Aw\, 

(5.10) W2 = BW2 

where A and B are linear operators (matrices). It is known that \i A = B then by 
linear diffusive coupling the two oscillators will synchronize and the synchronization 
is stable for large enough coupling strength. Now in the more general situation A 
and B can be different, either because the true systems are not made perfectly exact, 
or other reasons. But we still want to couple them so that they achieve identical 
synchronization. The analysis in the previous sections offer us a way to design the 
coupling form. 

Again use one way coupling and assume the form of coupling function on W2 ■ The 
coupled system becomes 

w\ = Awi, 

(5.11) W2 = Bw2 + <j{L{wi) — HW2) 

where _ff is a coupling matrix (known in advance) and L is the unknown coupling 
function of wi that needs to be determined. By Eq. (|5.4p and (|5.6p . we find that: 

(5.12) L{wi) ^ (^^^ + h\ wi. 

Here we see that only when A ~ B we can allow L{wi) = H{wi) and yield identical 
synchronization, otherwise the coupling form of x and y shall be different to compen- 
sate for the difference of the original dynamics on x and y respectively. The stability 
of the manifold is described by 

(5.13) i = [B-aH)^ 

following Eq. (|5.9p . In order for the synchronization manifold to be stable, we need 

(5.14) p{B~aH)<l, 

where p{Q) is the spectral radius of the matrix Q. 

Eq.s (|5.12[) . (|5.13[) and (|5.14p gives us a way to easily manipulate the coupling 
functions in the non-identical linear systems to obtain stable identical synchronization. 
With minor modifications we can obtain similar results for two-way coupling, and with 
a few more modifications as well as more analysis we shall be able to control more 
complicated systems. 
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6. Discussion and Conclusion. In this paper we have developed systematic 
methods for explicit construction of generalized synchronization manifolds in systems 
of coupled non-identical oscillators by means of the manifold equation, a PDE that 
must be satisfied by the manifold, and the associated variational equation describing 
its stability. 

Although the manifold equation gives necessary condition to determine the in- 
variant manifold, it is not sufhcient, mainly due to the fact that boundary conditions 
for this PDE are not easily attainable. To obtain a solution of the manifold equation 
without knowing explicitly the boundary conditions, we have proposed a time depen- 
dent PDE whose stably stationary solution is a solution of the manifold equation, 
and developed iteration schemes to solve it both symbolically and numerically. Sev- 
eral examples of constructing synchronization manifolds have been given for systems 
of nearly identical oscillators (considered as a perturbed version of identical oscil- 
lators, which correspond to complete synchronization) with unidirectional coupling, 
where we obtained the first order approximation of the perturbed manifold. A gen- 
eral technique for using the manifold equation to design the coupling for the purpose 
of controlling the form and stability of synchronization has also been discussed and 
illustrated using the simple case of linear systems. 

Our nunrcrical algorithms used to obtain the stationary solution of the time depen- 
dent PDE is analogous to the Euler method in a functional space, whose convergence 
is not necessarily global. Thus, finding a systematic method to choose an appropri- 
ate initial condition, as well as a detailed convergence analysis for the scheme, is an 
important problem that must be addressed in the future. As there exists a variety 
of different schemes to solve PDEs in the literature, exploring and comparing with 
other numerical schemes for solving the manifold equation is another important topic 
of future research. 
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Appendix. In the coupld Van der Pol system (|4.15p . in order to obtain the first 
order approximation of the synchronization manifold in e, we neglect higher order 
terms in e and write 

2^2 = a:^i + O(e^) « xi, 
(6.1) 2/2 = yi + eV{xi,yi) + 0{e^) « yi + eV{xuyi). 

Thus, using Eq. (|4.15[) . we have 

y2-yi=eil-xl)y,+Oie^). 
(6.2) 

We also have, for ^ <C 1, 

fdV. dV .\ ^,2, 
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Replacing xi, yi with x and y for convenience and equating the first-order e terms 
in the above two equations, we obtain the following equation: 

dV dV 

(6.4) y__^_ = (l_:,2)y 

This equation has a particular solution: 

(6.5) Vix,y)^x-^x^ 

Since the difference between x components in this case is already O(e^), we take 
the identical manifold X2 = xi as the first order approximation, while for the y 
component we use the form obtained above: 2/2 = 2/1 + e^(a^i, J/i)- 
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